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We discuss the possibility of forecasting earthquakes by means of (anti)neutrino tomography. Antineutrinos 
emitted from reactors are used as a probe. As the antineutrinos traverse through a region prone to earthquakes, 
observable variations in the matter effect on the antineutrino oscillation would provide a tomography of the 
vicinity of the region. In this preliminary work, we adopt a simplified model for the geometrical profile and 
matter density in a fault zone. We calculate the survival probability of electron antineutrinos for cases without 
and with an anomalous accumulation of electrons which can be considered as a clear signal of the coming 
earthquake, at the geological region with a fault zone, and find that the variation may reach as much as 3% for 
i'e emitted from a reactor. The case for a Ve beam from a neutrino factory is also investigated, and it is noted 
that, because of the typically high energy associated with such neutrinos, the oscillation length is too large and 
the resultant variation is not practically observable. Our conclusion is that with the present reactor facilities 
and detection techniques, it is still a difficult task to make an earthquake forecast using such a scheme, though 
it seems to be possible from a theoretical point of view while ignoring some uncertainties. However, with the 
development of the geology, especially the knowledge about the fault zone, and with the improvement of the 
detection techniques, etc., there is hope that a medium-term earthquake forecast would be feasible. 

PACS numbers: 14.60.Pq, 13.15.+g 91.30.Px 91.35.Pn 



Introduction. — Earthquakes and tsunamis are natural catas- 
trophes. They occur on our globe so frequently and cost thou- 
sands of lives and immense loss of property. We lack an effi- 
cient way to forecast earthquakes at present. We believe that 
our knowledge of modern physics and sophisticated detection 
techniques and facilities may help. It has been proposed that 
we could explore oil storage with neutrino beams Ll.-^], which 
also find applications in various other fields 10. S]. These cre- 
ative ideas motivate us to investigate the possibility of fore- 
casting earthquakes and tsunamis in terms of matter effects 
on (anti)neutrino oscillations. 

Various neutrino experiments are, or will soon be, in op- 
eration. With great effort, most neutrino parameters have 
been determined with certain accuracy. However, the complex 
phase Sep and the sign of Arri^i are still unknown. Moreover, 
so far we only have the upper bound of 613, which is man- 
ifestly smaller compared with the other two mixing angles. 
Precise measurement of 6*13 will be valuable for searches for 
CP violation in the lepton sector. Two reactor neutrino ex- 
periments, Daya Bay |6] in China and Double Chooz |7] in 
France, aiming to measure 6*13 are expected to reach very high 
precision at the percent level. 

Matter effects on the neutrino oscillation are crucial in some 
cases, for example, to explain the solar neutrino flux deficit. 
Wolfenstein formulated the matter effect on neutrino oscil- 
lation and then Mikheev and Smirnov [gj applied and devel- 
oped the theory to successfully solve the solar neutrino prob- 
lem. The newest experimental data confirm the MSW model. 

As a neutrino traverses through matter, its oscillation length 
is no longer the same as in vacuum. Electron neutrinos acquire 
an extra potential V{x) = ^/2G Fne{x) in the Hamiltonian 
due to their interaction with electrons, which exist in ordinary 
matter, where Gf is the Fermi coupling constant and n.e is 
the number density of electrons in matter. In a simplified ver- 
sion with only two generations, the expression of the survival 



probability in uniform matter is similar to the vacuum case, 
but the vacuum mixing angle and mass-squared difference 
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quantities in matter. That is 



mi I should be replaced by the corresponding 
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where 6m and Am|^ are the mixing angle and mass-squared 
difference in matter, respectively. and L are, respectively, 
the energy of the neutrino and the length of the baseline. For 
9 M we have 
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For antineutrinos, the minus sign in Eq. dJJl should be replaced 
by a plus sign, which implies that unlike the case of neutrino 
oscillation in the sun, no resonant stage exists, but an enhance- 
ment of the survival probability for antineutrinos, as discussed 
intlfl]. 

In analogy to the X-ray tomography, the possibility that 
(anti)neutrinos could be treated as a probe to detect the in- 
ner structure of the earth has been discussed by several au- 
thors. Two approaches | I 11] have been suggested for different 
(anti)neutrino energies. 

Our simplified model for an earthquake fault zone. — 
Earthquakes are caused by rapid slippage along faults and re- 
lease of energy stored in the seismic zone. Generally, a seis- 
mic zone does not consist of one single fault. The term " 
fault zone" implies that the zone is composed of several in- 
ner faults and undergoes tremendous deformation caused by 
severe shear and strain. 

On May 12 of 2008, an Ms= 8.0 earthquake occurred in 
Wenchuan, Sichuan Province, China. It was caused by the 
fracture of faults in the Longmenshan fault zone. This fault 
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FIG. 1 : The configuration of the trajectory of the (anti)neutrino beam 
which is produced in the source S, propagates through slabs I, II and 
III, and reaches the detector D. The baseline length L=Li+w+Lui- 
The geometry of the fault zone is described by three parameters, 
length /, depth d and width w. The matter density of each slab is 
also marked in the figure. 



zone is composed of three parallel faults. The Wenchuan- 
Maoxian fault is the west most one. According to geological 
studies, the Longmenshan fault zone is stressed by the west 
Bayan Har Block and Sichuan Basin 11211 . and a 300 km long 
and 50 km wide deformation zone was formed horizontally. 

The simplified geometric structure of the fault zone, which 
we are going to use in this work, is a cuboid with length I = 
300 km, width w = 50 km and depth d = IQ km, and with 
an upper edge parallel to the earth's surface. In general, a 
crust extends to a depth of 30 km U3i1 , and this cuboid is just 
located at its middle layer The main shock of the Wenchuan 
earthquake has a hypocenter depth of about 14 km, and this 
hypocenter is inside the cuboid. A schematic of the cuboid is 
shown in Fig. [T] 

The geological structure with a fault zone is not stable and 
sometimes varies actively. The drastic variations within the 
region of such geological structure can induce severe earth- 
quakes. Before an earthquake takes place, the variation in the 
matter density, or more exactly, the electron density of the 
fault zone is the most significant factor. The relevant effects 
are categorized into four aspects: 



Strain. Strain can change the matter density directly. 
During the Wenchuan earthquake, the Longmenshan 
fault zone is strained by the west and east geological 
blocks. The map of strain intensity drawn by a pro- 
fessional research group lll4ll shows that the eastern 
edge of the Qinghai-Xizang Plateau is close to the re- 
gion where strain intensity changes steeply. The Long- 
menshan fault zone undergoes the steepest change. The 
strain intensity on the west side is four times larger than 
that on the east side. 

Shear. Shear plays an important role in the slippage 
of the fault. It can either increase or decrease the matter 
density. From the maximum shear strain rate map of the 
Bayan Har Block and Sichuan Basin LlZl . one can ob- 
serve that the maximum shear strain rate of the east of 
the Longmenshan fault zone is 5xl0~^/a, whereas on 
the west side, Bayan Har Block, it is 15 ^ 50x 
The picture indicates that the shear contrast of both 
sides is 3 ^ 10. 

Anomalous electric field. The total electron content 



(TEC) anomaly, as one of the precursors of an earth- 
quake, has been observed from the TEC map using the 
GPS satellites. The authors of Ref. HI] point out 
that, for strong mid-latitudinal earthquakes, the seismo- 
ionospheric anomalies are manifested as a local TEC 
increase or decrease in the vicinity of the forthcom- 
ing earthquake epicenter, and this phenomenon appears 
several days prior to the main shock and the TEC mod- 
ification can be as large as 30 ^ 90 % of the normal 
value. The physical mechanism for such ionospheric 
anomalies is that an anomalous electric field is pro- 
duced near the epicenter. Prior to an earthquake, inten- 
sive gas discharges occurring at the crust in the earth- 
quake preparation zone lead to an anomalously strong 
vertical electric field Hill . The TEC disturbances are 
the vertical drift of the ionospheric plasma under the in- 
fluence of the anomalous electric field. On the other 
hand, if the anomalous vertical electric field definitely 
exerts influence on the earthquake preparation zone, it 
can also change the distribution of free electrons in the 
fault zone. In particular, the research into conductiv- 
ity characteristics indicates that a low resistance thrust 
fault zone in the middle layer of the Longmenshan fault 
zone was observed |17|] . 
• Rock dehydration. Under the influence of tremendous 
strain and shear, rock would undergo a dehydration pro- 
cess, which also contributes to the density shift in the 
fault zone li sll . 

The first two aspects concern the geometrical deforma- 
tion, which contributes an accumulation of mechanical en- 
ergy. When the shear exceeded the friction strength, the faults 
in the fault zone were unlocked and slipped suddenly jioll . 

Density shift and estimation of fault zone width. — As an- 
tineutrinos travel inside the medium, the matter potential 
V{x) is linearly proportional to the electron density. For the 
medium in the crust, it is reasonable to assume that ~ Np 
w Nn, i.e. the electron fraction number Ye = Ne/{Np + Nn) 
w 0.5, where Ne, Np and Nn correspond to electron, proton 
and neutron numbers in the crust layer, and write the electron 
density in the normal crust matter as 
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where po is the matter density of the normal medium and itin 
is the mass of nucleon. As discussed above, the electron den- 
sity is shifted in the fault zone as 



(4) 



As pre-earthquake activities in the fault zone are taken into 
account, p varies from po into an effective matter density p'^ 



Pf = VPo, 



(5) 



which means that in the fault zone the effective matter den- 
sity shortly before an earthquake, pj, is rj times larger than 
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the normal matter density pQ. In fact, because the existence 
of the anomalous electric field does not directly influence the 
real matter density, but only contributes to the electron density 
of the fault zone, the real matter density of the fault zone pf 
would be smaller than rjpo, i.e. 



2mNn{ = p'f ^ Wo, 
2mNn^ ^ Pf < vpo, 



(6) 



where is the number density of protons in the fault zone 
cuboid. The experimental observation indicates that there is a 
definite relation between p'y and pf, namely as increases, 
Pf increases too. The antineutrino oscillation determines p'p 
and with the data on p'p we would obtain unambiguous infor- 
mation about Pf which causes an earthquake. 

The study of geology focuses on a single fault or single 
fracture. Instead, we consider a zone that consists of several 
faults and has a relatively large width. The fault zone is treated 
as a slab with a uniform effective matter density p'p 

The width w of the fault zone can also be estimated from 
an empirical formula lfl6ll 
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FIG. 2: The differences (in percent) in survival probabilities of 
the reactor electron antineutrino between cases without and with an 
anomalous accumulation of electrons whose rj =3 and width w = 20, 
50, 80 and 100 km, respectively. 



for the long propagation ifioll . Then the evolution equation is 

ax 2h,, 



where Ms is the magnitude of the earthquake. For the 
Wenchuan earthquake, Ms = 8.0, the width w is about 80 
km from Eq. (|7]i, which is the fault zone width we use in our 
numerical computations. 

Configuration and result for the reactor antineutrinos. — 
The Daya Bay reactor located in Shenzhen is an example and 
we shall use its corresponding data as the inputs of our nu- 
merical computations. The fault zone cuboid in Wenchuan is 
about 1418 km to Daya Bay. As is shown in the configuration 
figure, the z/g first traverses the normal medium with a matter 
density pQ = 2.8 g/cm? and then enters the fault zone, which 
has a width w = 50 km. 

Since fault zones are not accessible for direct observation, 
all information concerning details of the zones is still coarse 
and sparse. Obviously, more concrete and precise informa- 
tion is needed very badly to forecast disasters. Along with the 
geologic study, theoretical work and new ideas are also neces- 
sary. Certainly, abundant geological data and development in 
seismology are helpful for the improvement of our fault zone 
picture and further research. In this work, we set r/ = 3 as a 
simple demonstration of our scheme. 

We suppose that the beam would be detected at the Qinghai- 
Xizang Plateau, which is 200 km away from the fault zone. 
This description of the total configuration of the trajectory is 
shown in Fig. [T] 

Under the approximation of ^ the vacuum oscilla- 
tion probability Pp^^p^ within the three-generation scenario 
can be reduced into the two-generation case. When taking the 
matter effect into account and considering the normal hierar- 
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'32 



the small mass-squared dif- 



ference Atoji dominates the oscillation as the effect induced 
by the larger mass-squared difference Am|j^ is averaged out 



where ijjakioi = e,p,T;k = 1, 2, 3) is the Pa 
evolution amplitude and M is given as ifioll 
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3,y — s'mOij, Cij = cos6ij and Acc = 2E^V{x). 
Under the approximation 6'i3 ^ 0, the amplitude i/)c3 is de- 
coupled from the evolution equations for ifjai and ilja2, and 



the effective survival probability 
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generation scenario just reduces into Pp 
pression in the two-generation case as 



(2) 



within the three- 
.r, which is the ex- 
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It is noted that the mixing angle 6*23 does not explicitly show 
up in the survival probability of D^- Since sin'^ ^13 < 

5 X 10-2 (99.73% CL) |20], with the assumption 9^^ - 0, the 
two-generation analysis involving two parameters Am^i and 
6*12 would be a good approximation for the realistic case. In 
fact, the mixing between and is sizable 121], therefore 
for the long-distance propagation, the oscillation between v^^ 
and Vr is almost complete, i.e. sin^ 26*23 ~ 1- Thus for the 
practical mixing between v^. and v^, really is the superposi- 
tions of and Vt- In the study of the solar neutrino flux miss- 
ing where only two generations are considered, the updated 
values of the input parameters are determined as sin^ 26*0 = 
0.86 and Ajti^ = 8x10"^ eV^. The e nerg y of the reactor 
electron-antineutrino of Ey = 3.6 MeV jioll is another input 
parameter for our numerical computation. The complete for- 
mulation involving three generations would make the whole 
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larger than one, the predicted probabihty difference undergoes 
a relatively larger change with the increase in rj. 

Neutrino factory. — Neutrino beams from a neutrino fac- 
tory i23[l can also be employed in the neutrino oscillation to- 
mography. In our configuration, we choose an electron neu- 
trino beam with energy = 1 GeV. Our calculation indi- 
cates that the difference in survival probabilities correspond- 
ing to the cases with and without the anomalous electron ac- 
cumulation is not obvious. Indeed, an observable difference 
of about 0.01% may be detected if the detector is 6000 km 
away from the source. However, such a length of baseline is 
much beyond the maximal distance in the crust layer and such 
a small difference may not be detected with the accuracy of 



FIG. 3: The dependence of the differences (in percent) in survival 
probabilities between cases without and with an anomalous accu- 
mulation of electrons on the rj value for reactor antineutrinos whose 
width w = 20, 50, 80 and 100 km, respectively. 



picture very complicated and does not change our qualitative 
conclusion. 

In this configuration, there are three slabs, I, II and III. Each 
of them is assumed to possess a uniform matter density and 
can be treated with an i?,-matrix 12211. 
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The elements of the i?-matrix are determined by the mass- 
squared difference Am|^, neutrino energy E,y and antineu- 
trino path length in the slab. In a uniform slab, the elements 
of the i?-matrices are 



Parameter 


Reactor 


Slab 


I 


II 


III 


L (km) 


1418 


w 


200 


P (ff/cm^) 


2.8 


2.8 X 3 


2.8 


Amli (10-'' eV^) 


8.02 


8.09 


8.02 


tan 26 m 


2.42 


2.30 


2.42 


Lfsc (km) 


111.18 


110.36 


111.18 


Parameter 


Neutrino Factory 


Slab 


I 


II 


III 


L (km) 


1418 


w 


200 


P (ff/cm^) 


2.8 


2.8 X 3 


2.8 


l^mii (lO'^eV^) 


1.46 


4.42 


1.46 


tan 2eM(10-^) 


-5.10 


-1.68 


-5.10 


L*L (lO^km) 


11.92 


3.92 


11.92 



TABLE I: The parameters of neutrino oscillation in matter for beams 
produced in a reactor and in a neutrino factory. For reactor electron 
antineutrinos Ei, = 3.6 MeV and for neutrino factory electron neutri- 
nos E^=7 GeV. 
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where x is the antineutrino path length in the slab, we cal- 
culate the survival probabilities without and with the anoma- 
lous accumulation of electrons at the fault zone. In Fig. |2] 
we present the differences in survival probabilities of the re- 
actor electron antineutrino between cases without and with 
an anomalous accumulation of electrons for 77 =3 and differ- 
ent widths. For oscillations in medium with and without the 
anomalous electron accumulation, the maximal difference in 
Pp^_>p_. is up to 2% when the width w = 100 km. 

Among all of the fault zone parameters, the width and the 
77 value play important roles in antineutrino oscillations. For 
reactor electron antineutrinos, the dependence of the differ- 
ences of survival probabilities on 77 is plotted in Fig. [3] with 
antineutrino energy E,y = 3.6 MeV for different width values. 
Fig. [3] indicates that for rj < 1, the effect is tiny, and for rj is 



the present neutrino facilities in our configuration. The reason 
is obvious that the accelerator neutrinos have much larger en- 
ergy than the reactor antineutrinos, thus the oscillation length 
igf J, for neutrinos from a neutrino factory is much longer than 
that for the reactor ones as shown in TableU 

Result from perturbation method. — In a low matter den- 
sity medium, V <C Ani^ /2E,y, loannisian and Smirnov 12411 
developed a perturbation method where e = 2VEi^/ Arri^ is 
taken as the expansion parameter for calculating the oscilla- 
tion probability. For the maximal matter density j]po in our 
case, one can determine 



2E^V 



0.029 < 1, 
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for El, = 3.6 MeV and 77 = 3. With the perturbation method, 
the exact adiabatic phase shift in matter is obtained as 
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FIG. 4: The difference (in percent) of survival probabilities between 
cases without and with an anomalous accumulation of electrons eval- 
uated in the perturbation theory Eq. l ll4t . 



and the plus and minus signs correspond to antineutrinos and 
neutrinos, respectively. By the formula for flavor-to-flavor 
transition, we obtain the survival probability of the electron 
antineutrino originating from the initial spot xq and ending at 
the final spot Xf. Up to the leading order of V{x), we have 



Pd^d 1 -sin2 26'o 
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In this scheme, we calculate the survival probability for re- 
actor electron antineutrinos. The result is shown in Fig. |4] 
The difference in survival probabilities for the cases with and 
without the anomalous electron accumulation can be as large 
as 3%, which is larger than the results obtained with the R- 
matrix method. 

Discussion. — We note that there is a sharp controversy 
about the rj value among geologists. Some think that this value 
should be rather small, whereas others advocate that, before a 
severe earthquake, the r/ might reach a sizable value. Taking 
the advice of the more conservative geologists, we adopt a rel- 
atively small value in this work. Indeed, we only investigate 
the physics and provide a possible way to forecast the earth- 
quake in terms of neutrino medium effects in this work, but we 
will leave the detailed study of the rj value to the geologists. 

In fact, the other neutrino sources, such as the cosmic 
muon, would contaminate our detection circumstance. How- 
ever, we may set some shielding and veto facilities to reduce 
the uncertainty caused by other neutrino sources. The detailed 
technical difficulties would be left for our experimental col- 
leagues. 

Foreseeing earthquakes in a certain region within ten years, 
two years and one or two months corresponds to long-, 
medium- and short-term earthquake forecast, while an im- 
pending earthquake forecast should be made a few days be- 
fore the outbreak. For the long-term and medium-term fore- 
cast, there is quite a long duration for the detector to collect 
sufficient data even though the flux is suppressed by the small 



solid angle spanned by the detector For short-term and im- 
pending forecasting, the electron density changes drastically 
as the earthquake is coming, and the overall effect results in a 
large shift in the oscillation probability in a short time inter- 
val, so that observing clear signals might be feasible if the flux 
from the source is sufficiently strong. As mentioned previ- 
ously, the geometrical deformation caused by strain and shear 
would last for a long period, say, hundreds of years in some 
cases and the deformation would eventually reach its maxi- 
mum before the outbreak. However, the strong anomalous 
electric field appears only several days, e.g. about a week 
before the earthquake, so that the week before the outbreak 
is crucial for detectors to collect data and reach a conclusion 
about the upcoming disaster 

Now let us make a rough numerical evaluation of the event 
rates based on the Daya Bay Monte-Carlo simulation. In the 
Daya Bay case, the detector is set at 1 km away from the neu- 
trino source, and the detection rate is about a few hundred per 
day (~400/day, as estimated). Supposing a shortened base- 
line in our case is 200km, the suppression factor for the flux 
would be 1/(200)^ ^ 3 x 10^^. Considering an exposure 
time of 100 days and assuming that the detection efficiency is 
100 times larger than that employed in Daya Bay experiment, 
the expected events would be at least 0(100), which is suit- 
able for a medium-term forecast. For example, in a period of 
one year, one can compare the events observed in the first half 
of the year with that in the second half. If these data show 
an explicit difference in survival probabilities during the two 
periods, they would imply a geological deformation that may 
result in an earthquake in the forthcoming period. The closer 
the earthquake is, the larger the difference. 

Indeed, in our present work, we consider idealized beam 
characteristics while ignoring many details. The uncertainty 
caused by other neutrino sources such as other reactors and 
geoneutrinos, is also ignored, which may be one of the main 
challenges in the implementation of this scheme. Thus we 
reach our conclusion that with the present reactor facilities 
and detection techniques, it is still a difficult task to make 
an earthquake forecast using such a scheme, though it seems 
to be possible from a theoretical point of view while ignor- 
ing some uncertainties. However, there is hope that with the 
great improvement in detection and innovation of facilities, 
the sensitivity and efficiency would be greatly improved and 
the size of the detector may be much larger with sufficient 
fund support, and with the development of the geology and 
seismology, the pictures of the fault zone will be clearer and 
its parameters more accurate, thus a medium-term earthquake 
forecast would be feasible. 

A scheme might remedy the shortcoming of remarkably 
losing neutrino flux due to the small solid angle. Namely, 
one may use a low energy beta beam 12511 as a source. Such 
a neutrino beam could be directly oriented to the fault zone. 
Without the 10^^ suppression factor, the requirement for the 
size and sensitivity of the detectors would be greatly allevi- 
ated. However, in order to obtain intense and coUimated neu- 
trino beams from the traditional source for the beta beam, e.g.. 
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f He++^-^Li+++ +e~ +Ve iH, the high intensity of the 
jHe^"'' source and a large Lorentz boost factor are necessary. 
We suggest alternative schemes. One is that in the process 
e~ + p — > n + j/c we can employ the very coUimated electron 
beam to bombard the protons at rest. Another is that we can 
use the lowly accelerated pion beam for tt^ — > + De{ve)- 
Even though in the second proposal, the decay rate is sup- 
pressed by the helicity rule, the branching ratio is still of the 
order of 10^^ [26], so with a large number of pions, it is 
still plausible. Moreover, a large Lorentz boost factor |p|/TOe 
which determines the spreading solid angle, is easy to obtain 
and the energy for the pion accelerator needs to be just at 
©(GeV). Thus the size of such facilities can be small and 
the corresponding technique would be relatively simple. 
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very enlightening discussions on the subject with Prof. Chung 
Ngoc Leung. In fact, this work was motivated by the discus- 
sion when he visited the Nankai campus. This work is sup- 
ported by the National Natural Science Foundation of China 
(NNSFC) and the Special Grant of the Education Ministry of 
China for the Ph.D programs. 
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